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Abstract 

The mean first exit time and escape probability are utilized to 
quantify dynamical behaviors of stochastic differential equations with 
non-Gaussian a— stable type Levy motions. Both deterministic quanti- 
ties are characterized by differential- integral equations (i.e., differential 
equations with nonlocal terms) but with different exterior conditions. 
The non-Gaussianity of noises manifests as nonlocality at the level of 
mean exit time and escape probability. An objective of this paper is to 
make mean exit time and escape probability as efficient computational 
tools, to the applied probability community, for quantifying stochastic 
dynamics. An accurate numerical scheme is developed and validated 
for computing the mean exit time and escape probability. Asymptotic 
solution for the mean exit time is given when the pure jump measure 
in the Levy motion is small. 

From both the analytical and numerical results, it is observed that 
the mean exit time depends strongly on the domain size and the value 
of a in the a— stable Levy jump measure. The mean exit time can 
measure which of the two competing factors in a-stable Levy motion, 
i.e. the jump frequency or the jump size, is dominant in helping a 
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process exit a bounded domain. The escape probability is shown to 
vary with the underlying vector field (i.e., drift). The mean exit time 
and escape probability could become discontinuous at the boundary 
of the domain, when the process is subject to certain deterministic 
potential and the value of a is in (0, 1). 

Key Words: Stochastic dynamical systems; non-Gaussian Levy 
motion; Levy jump measure; First exit time; double-well system 

Mathematics Subject Classifications (2000): 60H15, 60F10, 
60G17 



1 Motivation 

Random fluctuations in complex systems in engineering and science are often 
non-Gaussian j28|. [9l [10]. For instance, it has been argued that diffusion 
by geophysical turbulence [26] corresponds, loosely speaking, to a series of 
"pauses" , when the particle is trapped by a coherent structure, and "flights" 
or "jumps" or other extreme events, when the particle moves in the jet flow. 
Paleoclimatic data [TT] also indicate such irregular processes. 

Levy motions are thought to be appropriate models for non-Gaussian 
processes with jumps [23j . Recall that a Levy motion L{t), or Lt, is a 
stochastic process with stationary and independent increments. That is, for 
any s, t with < s < t, the distribution of Lt — Lg only depends oiit — s, and 
for any < to < < • • • < ^n, -^^t, — ^u^i, ^ = !> ' ' ' j't-; are independent. 
Without loss of generality, we may assume that the sample paths of Lt are 
almost surely right continuous with left limits. 

This generalizes the Brownian motion B{t), which satisfies all these three 
conditions. But additionally, (i) almost all sample paths of the Brownian 
motion are continuous in time in the usual sense and (ii) Brownian motion's 
increments are Gaussian distributed. 

Stochastic differential equations (SDEs) with non-Gaussian Levy noises 
have attracted much attention recently |2l[23]. To be specific, let us consider 
the following scalar SDE with a non-Gaussian Levy motion 

dXt = f{Xt)dt + dLt, Xo = x, (1) 

where / is a vector field (or drift), and Lt is a scalar Levy motion defined 
in a probability space (ri, P). 

We study the first exit problem for the solution process Xt from bounded 
domains. The exit phenomenon, i.e., escaping from a bounded domain in 
state space, is an impact of randomness on the evolution of such dynamical 
systems. Two concepts are applied to quantify the exit phenomenon: mean 
exit time and escape probability. 

We define the first exit time from the spatial domain D as follows: 

t{uj) := mf{t > 0,Xt{uj,x) ^ D}, 
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and the mean exit time is denoted as u{x) := Er. The hkehhood of a particle 
(or a solution path Xf), starting at a point x, first escapes a domain D and 
lands in a subset E of V (the complement of D) is called escape probability 
and is denoted as Pe{x). 

The existing work on mean exit time gives asymptotic estimate for u{x) 
when the noise intensity is sufficiently small, i.e., the noise term in ([1]) is 
e dLt with < e ^ 1. See, for example, Imkeller and Pavlyukevich [151 116j . 
and Yang and Duan [29] , 

In the present paper, however, we numerically investigate mean exit time 
and escape probability for arbitrary noise intensity. The mean exit time 
u{x) and escape probability Pe{x) for a dynamical system, driven by a non- 
Gaussian, discontinuous (with jumps) Levy motion, are described by two 
similar differential-integral equations but with different exterior conditions. 
The non-Gaussianity of the noise manifests as nonlocality at the level of the 
mean exit time and escape probability. We consider a numerical approach 
for solving these differential-integral (non-local) equations. A computational 
analysis is conducted to investigate the relative importance of jump measure, 
diffusion coefficient and non-Gaussianity in affecting mean exit time and 
escape probability. 

Our goal is to make mean exit time and escape probability as efficient 
computational tools, to the applied probability community, for quantifying 
stochastic dynamics. 

This paper is organized as follows. In section 2, we recall the generators 
for Levy motions. In section 3, we consider SDEs driven by a combination of 
Brownian motion and a symmetric a— stable process. Numerical approaches 
and simulation results are presented in section 4 and 5, respectively. Finally, 
the results are summarized in section 6. 

2 Levy motion 

A scalar Levy motion is characterized by a linear coefficient 0, a diffusion 
parameter d > and a non-negative Borel measure defined on (M, ;B(M)) 
and concentrated on M \ {0}, which satisfies 



This measure u is the so called the Levy jump measure of the Levy motion 
L{t). We also call {9,d,v) the generating triplet. 

Let Lt be a Levy process with the generating triplet {6, d, v). It is known 
that a scalar Levy motion is completely determined by the Levy-Khintchine 





or equivalently 




(3) 
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formula (See [21 [231 [22]). This says that for any one-dimensional Levy pro- 
cess Lt, there exists a 9 £ R, d > and a measure such that 

Ee'^'^' = exp{ieXt -dt^+t [ (e*^^ - 1 - iXyI{\y\^,})u{dy)}, (4) 

^ Jr\{o} 

where Is is the indicator function of the set S, i.e., it takes value 1 on this 
set and takes zero value otherwise: 



1, ifyeS; 
0, iiy^S. 



The generator A of the process Lt is defined as Aip = limj|o '^"^"'^ where 
Pt(p{x) = Ex(p{Lt) and (p is any function belonging to the domain of the 
operator A. Recall that the space C^(K) of functions with bounded 
derivatives up to order 2 is contained in the domain of A, and that for every 
if G C^{m) (See [21 [22]) 

Aip{x) = eip'{x) + -ip"{x) + / [ip{x + y)- '^{x) - I{\y\<i} yp'ix)] v{dy). 
^ JlR\{0} 

(5) 

Moreover, the generator for the process Xt in ([T|) is then 

A^ = f(x)ip\x) + eip'{x) + ^^"{x) 

+ / [^{x + y)-^{x)-I{\y\<i}yip'{x)]u{dy). (6) 
Jr\{o} 

For a G (0,2], a symmetric a-stable process is a Levy process Lt such 
that 

A symmetric 2-stable process is simply a Brownian motion. When a £ (0, 2), 
the generating triplet of the symmetric a-stable process Lt is (0, 0, f^), where 

Ua{dx) = C„|2;|-(^+°) dx 

Q, P(l+^* ) 

with Cn, given by the formula C„ = — ; p= — — ^—r— . For more information 

"B .r a 2i-°0Fr(l-f) 

see 



[Il[8l[2l]. 



3 Mean exit time and escape probability 

Now we consider the SDE ([1]) with a Levy motion Lf that has the generating 
triplet (0,(i, efc), i.e., zero linear coefficient, diffusion coefficient d > and 
Levy measure ei'aidu), with < a < 2. This Levy motion is the independent 
sum of a Brownian motion and a symmetric a-stable process. Here e is a 
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non-negative parameter and it does not have to be sufficiently small. Strictly 
speaking, this is not a a— stable Levy motion (because the diffusion d may 
be nonzero), but a Levy motion whose jump measure is the same as that of 
a a— stable Levy motion. 

We first consider the mean exit time, u{x) > 0, for an orbit starting at 
X, from a bounded interval D. By the Dynkin formula [21 [23] for Markov 
processes, as in [HI [251 EI], we obtain that u{x) satisfies the following 
differential-integral equation: 

Au{x) = -1, xe D (7) 
n = 0, X £ D", 

where the generator A is 

Au = f{x)u (x) + 



+e / [u{x + y)- u{x) - Iiiyi^n yu'{x)] i^a{dy), (8) 
Jr\{o} 

and D'^ = M\D is the complement set of D. 

We also consider the escape probability of a particle whose motion is 
described by the SDE ([T]). The likelihood of a particle, starting at a point x, 
first escapes a domain D and lands in a subset E of (the complement of 
D) is called escape probability. This escape probability, denoted by Pe{x), 
satisfies [18\ [27] the following equation 

APe{x) = 0, x£D, (9) 
PeIxge = 1, Pe\x&D''\E = 0> 

where A is the generator defined in ([8|). 

In the present paper, we only consider scalar SDEs. For SDEs in higher 
dimensions, both mean exit time and escape probability will satisfy partial 
differential- integral equations, and our approaches generally apply. 



4 Numerical schemes 



i{\y\<s}iy)y 



Noting the principal value of the integral / — dy always vanishes 

"/ M 1^1 

for any J > 0, we will choose the value of 6 in Eq. ([7|) differently according 
to the value of x. Eq. ((Tj) becomes 

- u {x)+f{x) u {x)+eCa / 1 M , ' dy = -1, 

2 Jr\{o} |yr+" 

(10) 

for X G (a, 6); and u{x) = for x ^ (a, b). 
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Numerical approaches for the mean exit time and escape probabihty in 
the SDEs with Brownian motions were considered in [6l [7], among others. 
In the following, we describe the numerical algorithms for the special case of 
{a,b) = (—1,1) for clarity of the presentation. The corresponding schemes 
for the general case can be extended easily. Because u vanishes outside 
(-1, 1), Eq. ([ini) can be simplified by writing = J'^"" + J^~^^ + J^^, 



d 



u"{x) + f{x)u'{x) 



a 



1 



+ 



1 



(l + x)" (1-2;) 



u{x) 



l-x 



u{x + y) - u{x) - I{\y\^s}yu'{x) 



l-x 



\l+a 



dy 



■1, (11) 



for x G ( — 1, 1); and u{x) = for x ^ ( — 1, 1). 

Noting u is not smooth at the boundary points x ■ 
ensure the integrand is smooth, we rewrite Eq. (jlip as 



-1, 1, in order to 



d eC 
-u"{x) + f{x)u'{x) 



+eCa 



a 



1 



+ 



1 



-l-x 

for X > 0, and 



\y 



l+a 



(1 + x)" (l-x) 
u{x + y) — u{x) — yu'{x) 



l+x 



\y 



l+a 



u{x) 
dy 



^u"{x) + /(x)n'(x) - ^ 
2 a 



1 



+ 



1 



+eCa 



l-x 



l+x 



u{x + y) — u{x) 



\y 



l+a 



dy + eCa 



(l + x)" (l-x 
u{x + y) — u{x) — yu'{x) 

-l-x \y 



u[x) 



l+a 



dy 



for X < 0. We have chosen 5 = min{| — 1 — x|, |1 — x|}. 

Let's divide the interval [—2, 2] into 4 J sub-intervals and define Xj = jh 
for — 2J < j < 2J integer, where h = 1/J. We denote the numerical solution 
of u at Xj by Uj. Let's discretize the integral-differential equation (fT2]l using 
central difference for derivatives and "punched-hole" trapezoidal rule 



d Uj^i 
2 



2Ui + 



/l2 

-J+j 



U. 



U. 



2h 



a 



1 



+ 



1 



(l + x,)" (l-x,) 



+ eCc^h 

k=-j-j 



+ eCah 



\Xk 



l+a 



J-j 

E 



„ Uj+k - Uj - {Uj+i - Uj^i)xk/2h 



\Xk 



l+a 



(14) 



where j = 0,l,2,---,J — 1. The modified summation symbol ^ " means 
that the quantities corresponding to the two end summation indices are 



multiplied by 1/2. 



dUj^i-2Uj + Uj+i Uj+i 



u. 



Xk 



l+a 



2 ^2 

k=J+j 

where j - 

values of Uj vanish if the index 



2h 



a 



1 



+ 



1 



+ (1-x,) 



J+j 



„ Uj+k - Uj - {Uj+i - Uj^i)xk/2h 



\Xk 



l+a 



(15) 



J + I,-- - ,— 2,— 1. The boundary conditions require that the 
> J- 

The truncation errors of the central differencing schemes for deriva- 
tives in ()14p and (jlSp are of 2nd-order 0{h?). From the error analysis of 
Navot (1961) pO], the leading-order error of the quadrature rule is —C,{a — 
l)u" {x)h?~°^ + O^h?), where Q is the Riemann zeta function. Thus, the fol- 
lowing scheme have 2nd-order accuracy for any < a < 2, j = 0, 1, 2, • • • , J— 
1 



U. 



i-1 



h? 
-J+j 
+ eCch ^ 



2h 



a 



1 



+ 



1 



(l+x,)° (l-x,) 



U, 



-J-j 



\Xk 



+ eCa^h Y 



„ Uj+k - Uj - {Uj+i - Uj-i)xk/2h 



-J+3,k^0 



\Xk 



l+a 



where ^h = — £C'aC(o — 1)^^ Similarly, for j = — J + 1, ■ 



u. 



k=j+j 



2h^ 

^11 Uj+k - Uj 



2h 



a 



+ 



(16) 

-2,-1, 
1 



{i + xjY {i-xjY 



u 



\Xk 



+ eCah Y 



„ Uj+k - Uj - {Uj+i - Uj^i)xk/2h 



-J-j,k^o 



\Xk 



l+a 



(17) 



where j = -J +!,■■■ , -2, -1, 0, 1, 2, • • • , J - 1. Uj = if \j\ > J. 

We solve the linear system (|16m7p by direct LU factorization or the 
Krylov subspace iterative method GMRES. 

We find that the desingularizing term {I^y^^s}yu' (x)) does not have any 
effect on the numerical results, regardless whether we use LU or GMRES 
for solving the linear system. In this case, we can discretize the following 
equation instead of PT|) 



-u"{x) + f{x)u'{x) 





1 


a 


Xl + x) 




nl—X 


+eC 


3! / 

J-l-X 



+ 



1 



u{x + y) — u{x) 



\y 



l+a 



u{x) 
dy 



(18) 
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where the integral in the equation is taken as Cauchy principal value integral. 
Consequently, instead of (jl6p and (jl7p , we have only one discretized equation 
for any < a < 2 and j = - J + 1, - ■ ■ , -2, -1, 0, 1, 2, • • • , J - 1 



C, 



Uj-i - 2Uj + Uj+i 



sCf^ Uj 
a 



/l2 



1 



U. 



U. 



+ 



+ (1-x,) 



2h 

+ eCah 



E 



-J-j,k^O 



\l+a 



(19) 



With minor changes, we also have the scheme for simulating escape prob- 
ability characterized by the equation Q. 



5 Numerical results 
5.1 Verification 

5.1.1 Comparing w^ith analytical solutions 

In order to verify that the numerical integration scheme for treating the 
improper integral in ([7j) is implemented correctly, we compute the left-hand 
side(LHS) of ([TD by substituting u{x) = I - x^, d = 0, f {x) = e = 1 and 
(a, 6) = (-1,1) 



LHS 



(1+x 



— Ca 
+ (1 

-Ca (4 + 2a;lni 



l-a / 1-x 
a 



2x _|_ 1+x 
1-a ~^ 2-a 



,1— a / 1+x I 2x I 1—x 



l-x 

X 



ifa^l; (20) 
if a = 1. 



Figure [T] shows the differences (the errors) between the numerical and the 
analytical values of LHS of ([7]) at the fixed value x = —0.5 for different reso- 
lutions J = 10, 20, 40, 80, 160 and 320. We plot log^g (error) against logio(^) 
for a = 0.5, 1, 1.5, where h = 1/J and error is the difference between the 
numerical and the analytical values of LHS. Clearly, the numerical results 
show that the error of computing LHS decays as 0{h?'). The second-order 
accuracy is expected from the error analysis of the numerical integration 
method ()19p . For a fixed resolution /i, the error increases as a increases due 
to the fact that LHS in (p0|) becomes more singular at x = — 1 and 1 as a 
increases. 

Next, we compare the numerical solution with the analytical solution for 
the mean exit time 1131 



u{x) 



2'^T{1 + a/2)T{l/2 + a/2) 



(21) 
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Figure 1: The error of the numerical values of the left-hand side of Eq. ^ 
compared with the analytical expression in for u{x) = 1 - x^, d = 0, 
f[x) = 0, e = 1 and (a, h) = (—1, 1). The results are computed at x = —0.5 
for different values of a = 0.5 (marked by *'s), 1 (the x's) and 1.5 (the o's) 
and different resolutions J = 10, 20, 40, 80, 160 and 320. Also shown is an 
illustrating solid line with slope equal to —2. 

in the special case of ([7]) in which ti = 0, f{x) = 0, e = 1 and (o, b) = {—b, b) 
with 6 > 0. Figure [2] shows the numerical solutions (the dashed lines) 
obtained by solving the discretized equations (jl6p and (jl7p or (|19p with the 
fixed resolution J = 80, (a, b) = (—1, 1) and different values of a = 0.5, 1, 1.5, 
while the corresponding analytical solutions are shown with the solid lines. 
The comparison shows that the numerical solutions are very accurate as 
one can hardly distinguish the numerical solution from the corresponding 
analytical one. Note that, in this case of {a,b) = (—1, 1), for a fixed value 
of the starting point x, the mean exit time u{x) decreases when a increases 
in the interval (0,2]. Later, we will see the dependence of the mean exit 
time on a is much more complicated when the size of the interval 6 — a is 
increased. 

Figure [3] gives the error in the numerical solution to the discretized equa- 
tions (fn|) and (fT5|) derived by using " punched- hole" trapezoidal rule. Here, 
we compute error = \u{—0.5) — C/(— 0.5)| by comparing the mean exit time 
at x = —0.5 for different resolutions and values of a, where u and U denotes 
the analytical and the numerical solution respectively. For a fixed resolution, 
the numerical error has similar sizes for a = 0.5 and a = 1 but it is much 
larger in the case of a = 1.5. The analysis shows that the rate of decay in 
the error as the resolution increases is 0(/i^~"). Our numerical results in 
Fig. [3] show slower decaying rates than those in the theory for a = 0.5 and 
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Figure 2: The mean exit time u{x) for the special case of the symmetric a- 
stable process d = 0, f{x) = 0, e = 1 and (a, 6) = (—1, 1). The dashed lines 
are the numerical solutions for different values of a = 0.5, 1, 1.5 obtained by 
solving the discretized equations (I19p with the resolution J = 80, while the 
solid lines represent the corresponding analytical solutions (|2ip including 
a = 2. 
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Figure 3: The error of the numerical solution of the mean-exit time u{x) from 
the discretized equations and ([15]) using the " punched- hole" trapezoidal 
rule for d = 0, f{x) = 0, e = 1 and {a,b) = (—1,1). The results are 
computed at x = —0.5 for different values of a = 0.5 (marked by *'s), 1 
(the x's) and 1.5 (the o's) and different resolutions J = 10,20,40,80,160 
and 320. Also shown are two illustrating solid lines with slope equal to — 1 
and —0.5 respectively. 



11 




Figure 4: The error of the numerical solution of the mean-exit time u(x) 
from the discretized equations (jl6p and (jl7p or (jl9p with the correction term 
for d = 0, /(x) = 0, e = 1 and {a,b) = (—1, 1). The results are computed 
at X = —0.5 for different values of a = 0.5 (marked by *'s), 1 (the x's) and 
1.5 (the o's) and different resolutions J = 10,20,40,80,160 and 320. Also 
shown is an illustrating solid line with slope equal to —1. 

1, which is due to the non-smoothness of the solution at x = —1, 1. 

Figure S] is the same as Fig. [3] except that the numerical results are 
obtained from the discretized equations ([T6|) and pT|) or ([T9|) with the cor- 
rection term that removes the leading-order quadrature error. Although 
the numerical analysis predicts the decaying rate of the numerical error is 
0(/i^), the numerical results shown in Fig. [4] indicate the rate of decay is 
only 0{h), because the analytical solution u{x) given in (|2ip has infinite 
derivatives at x = — 1 and 1. Note that we have demonstrated in Fig. [T] that 
the convergence order would be 2 if the solution u were smooth on the whole 
closed interval [—1,1]. Though the convergence order is only 1, it become 
independent of a after we add the correction term and the numerical error 
is two orders of magnitude smaller than that without the correction term 
when a = 1.5. 

5.1.2 Comparing with asymptotic solutions 

Since we do not have analytical solutions in closed form for the mean first 
exit time u to Eq. ^ in the general case of d 7^ 0, / 7^ 0, we calculate the 
asymptotic solutions for small values of e, i.e., small pure jump measure in 
the Levy motion. Then, we test our numerical schemes by comparing the 
results with those from the corresponding asymptotic solutions. 
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We look for solution to Eq. ([7]) in the form of expansion 

u{x) = uo{x) + eui{x) + e'^U2{x) H . (22) 

By setting e = in Eq. ([7|), we obtain the equation for uq 

^4(x) + fix)u'oix) = -l. (23) 

The solution is given by 

r r f^y^'^y dz + Ai 

uo{x)= v{y)dy + A2, where v{x) = 2 , . ,(24) 

where Ai and A2 are integration constants that can be determined by the 
boundary conditions. Substituting the expansion ([22]) into d?]) and discard- 
ing the terms with second or higher powers of e, we obtain the equation for 

Ul 

-u^{x)+f{x)ui{x)+Ca ^ dy = 0. 

(25) 

^. ^ f + y) - M^) - i{\y\<i}iy)yuo{x) 

Denotmg g{x) := Ca / ^ |1+q — '^^ have 



(26) 

where A3 and A4 are integration constants that can be determined by the 
boundary conditions. 

Let us consider the special case /(x) = 0, d = 1 and (a, 6) = (—1,1). 
According to the general solution ()24p and ()26p . the zeroth and first-order 
solutions are 



uo{x) = l-x"^, (27) 



ui{x) 



[(2 -«)(! + 



a(l-a)(2-o)(3-a)(4-a) 

+(4 - q)(1 + x)3-°(l _ x) + (4 - a)(l + x)(l - x)^-" 

+(2 - a)(l - x)^-" - 24-"(2 - a)] , for a G (0, 1) U (1, 2), 

^ [2x2-2-41n2+(2 + x)(l-x)2ln(l-x) 

+ (2-x)(l+x)2ln(l+x)] , fora = l. 



To further verify our numerical methods, we compare the numerical so- 
lution to d?]) with the asymptotic solution ([27|) uq + £Ui for the case e = 0.1, 
/ = 0, d = 1 and (a, 6) = (—1, 1), as shown in Fig. \5\ The plots show that 
the two solutions are very close for a = 0.5 (Fig. [5]^a)), a = 1 (Fig. [5]^b)) 
and a = 1.5 (Fig. [5]^c)). Furthermore, Fig. [5]^d) shows that the difference 
between the two solutions is proportional to e^, which is expected as the 
asymptotic solution uq + eui given by (j271) is only accurate up to 0(e). 
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Figure 5: Comparison between the numerical solution to ([7]) and the asymp- 
totic solution (j27|) for small e with / = 0, d = 1 and {a,b) = (—1, 1). (a) 
e = 0.1 and a = 0.5. The numerical solution is displayed with dashed line 
while the asymptotic solution uq + eui is shown with solid black line, (b) 
Same as (a) except a = 1. (c) Same as (a) except a = 1.5. (d) The dif- 
ference between the numerical solution U (0) and the asymptotic solution 
uo{0) + eui{0) is plotted against for a = 1.5 and the fixed resolution 
J = 100. 
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Figure 6: Dependence of the mean exit time u{x) on the size of domain 
D for pure a-stable Levy motion, i.e., d = 0, f{x) = 0, e = 1. (a) The 
mean first exit time u{x) from the analytical solution (12 ip for the domain 
D = (-1.5,1.5) and a = 0.5,1,1.5,2. (b) Same as (a) except D = (-4,4). 

5.2 Dependence of the mean exit time on the size of domain 
5.2.1 Pure jump: d = 0, f = 

It is well-known that the a-stable Levy process has larger jumps with lower 
jump frequencies for small values oi a (0 < a < 1) while it has smaller 
jumps with higher jump probabilities for values of a closer to 2. Given 
a bounded domain D, it would be interesting to know, for symmetric a- 
stable Levy motion (d = and / = 0), the exit times out of domain D are 
shorter for small values of a or large values of a. The answer, given by the 
analytic solution (12ip . depends on the size of domain D. To illustrate the 
dependence, we compare the solutions for four representative values of a, 
i.e., 0.5, 1, 1.5, and 2. For a small-size domain D = (—6, b) with < 6 < 1.25, 
it is easier to leave the domain for larger values of a and any starting point 
X, such as the case 6 = 1 shown in Fig. [2j In contrast, for large domains D 
with 6 > 3, the exit times are shorter when a is smaller except for starting 
points near the boundary such as the case b = 4 shown in Fig. [6]^b). For 
median-sized domains D = {—b,b) with 1.25 < 6 < 3, whether the mean 
exit time is shorter for smaller or larger values of a depends on both the 
position of the starting point x and the size of symmetric domain 26, such 
as the case b = 1.5 shown in Fig. [6)^a). 

The results on the mean exit time show that one has to consider both the 
domain size and the value of a when deciding which of the two competing 
factors in a-stable Levy motion, the jump frequency or the jump size, is 
dominant. The small jumps with high frequency, corresponding to Levy 
motion with a closer to 2, make it easier to exit the small domains. On 
the other hand, it is easier to exit large domains for a-stable Levy motion 
with a closer to 0, which has the characteristics of the large jumps with 
low frequency. Another observation from the results is that, for smaller 
values of a (0 < a < 1), the mean exit time profiles becomes flatter away 
from the boundary, such as the graphs for a = 0.5 in Fig. [2] and Fig. [6]^b). 
This implies that the jump sizes of the processes are usually larger than the 
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Figure 7: Dependence of the mean exit time u{x) on the size of domain 
D for the case of Ornstein-Uhlenbeck potential {f{x) = —x) with both 
Gaussian (d = 1) and non-Gaussian (e = 1) noises, (a) The mean first exit 
time u{x) for D = (—1,1) and a = 0.5,1,1.5,2. (b) Same as (a) except 
D = (—1.5,1.5). (c) Same as (a) except D = (—2,2). (d) Same as (a) 
except D = (-4,4). 

domain sizes, thus the mean exit times have small variations for different 
starting positions. 

5.2.2 Ornstein-Uhlenbeck(O-U) potential: f{x) = —x 

In the deterministic case dXt = —Xt dt, the origin is the sole stable point 
and the particle is driven toward the origin with the velocity proportional 
to its distance to the unique stable point. When a particle is subject to the 
Levy motion ([T|) defined in the beginning of Sec. [3] the mean exit time out 
of a bounded domain becomes finite. We emphasize that, in this paper, we 
consider the Levy motion defined by the generator in ([8|) where we vary the 
diffusion coefficient d and the parameter e independently. Figure[7]illustrates 
the dependence of the exit time u on the size of the domain for the case 
f{x) = —X, d = 1, e = 1 and the four typical values of a = 0.5, 1, 1.5,2. 
Figure El^a) shows that, for small domains D = (—6, 5) such as 6 = 1 and 
from any starting point x in Z), the mean exit times are shorter for larger 
values of a, in agreement with the corresponding result in the absence of 
the driving force / = and the Gaussian noise d = Q shown in Fig. [2j On 
the other hand, for large domains such as D = (—6, h) with 6 = 4, FiglTl^d) 
demonstrates that the mean exit times are longer for larger values of a. 
Again, the behavior agrees in general with the results in the previous pure 
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Figure 8: Dependence of the mean exit time u{x) on the size of domain D 
for the case of Ornstein-Uhlenbeck potential (/(x) = —x) with pure non- 
Gaussian noises d = and e = 1. (a) The mean first exit time u{x) for 
D = (—1, 1) and a = 0.5, 1, 1.5, 2. (b) Same as (a) except D = (—1.5, 1.5). 
(c) Same as (a) except D = (—2,2). (d) Same as (a) except D = (—4,4). 

jump case with f = d = and e = 1 shown in Fig. E^b). Note that, for 
the starting points near the boundaries, the relations between the mean exit 
times and the values of a are different in these two cases: the mean exit 
times increases as a is raised for the case of nonzero Gaussian noise d = 1 
and the nonzero driving force f{x) = —x while the pure jump Levy motion 
(/ = and d = 0) has the opposite behavior. The dependence of the mean 
exit times on the value of a is mixed for median-sized domains, such as 
D = (—1.5, 1.5) and D = (—2,2) shown in FigiT^b) and (c) respectively. 

Figure [8] shows the mean first exit times when Gaussian noise is removed 
while keeping other factors the same, i.e., f{x) = —x, e = 1 but d = 0. The 
dependence on the size of the domain is similar to the previous case with 
Gaussian noise d = 1. However, as shown in Fig. [HI we find that, in the 
presence of the 0-U potential {f{x) = —x) and without Gaussian noise 
(d = 0), the mean first exit time u{x) not only has a flat proflle also is 
discontinuous at the boundaries x = ±6 for a = 0.5. We find that it is 
also true for other values of a in (0, 1). A possible explanation for this is 
as follows: When a G (0,1), the original first order differential operator 
plays the dominant role, while when a G (1,2), the integral operator plays 
the dominant role. To obtain the discontinuous numerical solutions in these 
cases, we replace the central differencing scheme for the term f[x)u'{x) in 
()19p with a second-order one-sided difference for the first and last interior 
grid point. 
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Figure 9: Effect of Gaussian noise on the mean exit times of the Levy motion 
L". (a) The mean exit times for f{x) = —x,e = l,a = 0.5, D = (—1,1) and 
d = 0, 0.1, 0.5, 1. (b) Same as (a) except a = 1.5. 



5.3 Effect of the noises 

Having discussed the dependence of the mean exit times on the value of 
a in the Levy motion Lf defined in the beginning of Sec. O we examine 
the effect of Gaussian noise on the profile of the mean exit time u{x) as a 
function of the location of the starting point x. We vary the values of the 
diffusion coefficient d and the parameter e in Levy measure independently. 
Thus, Lf in the SDE, as defined by the generator in ([8|), is different than the 
traditional Levy motion where the diffusion coefficient d and the coefficient 
e in Levy measure are changed in tandem. 

Consider the Levy motion driven by 0-U potential f{x) = —x and the 
symmetric domain D = (—1, 1). Figure [9] shows the numerical results of 
the mean exit times with e = 1 for different amount of Gaussian noises 
d = 0,0.1,0.5,1. For small values of a (0 < a < 1), such as a = 0.5 as 
shown in Fig. [91(a), the mean exit time shapes as function of the starting 
point x change dramatically as the amount of Gaussian noises d increases. 
For small amount of Gaussian noises, the mean exit time profile is flat in the 
middle of the domain and drops to zero quickly near the boundary points; 
for large amount of Gaussian noises, the mean exit time profile become more 
parabla-like as shown in the graph for d = 1. 

It is worth pointing out that the mean first exit time is discontinuous 
at the boundary x = ±1 in the case of pure non-Gaussian noise d = and 
a € (0,1), i.e., the limits limx^±iu{x) are nonzero while it(±l) = 0. As 
mentioned in previous section, we have to use an one-sided difference scheme 
near the boundary to avoid numerically differentiating across discontinuities. 
From our numerical simulations for other values of a and domain sizes (not 
shown here), we find that the mean exit time u{x) driven by 0-U potential 
with "pure" a-stable jump only and < a < 1 would be discontinuous at the 
boundary of the domain. Recall that, in the absence of deterministic driving 
force (/ = 0) and Gaussian noise {d = 0), the mean exit time profile given 
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Figure 10: Effect of non-Gaussian noise on the mean exit times of the Levy 
motion defined by the generator in ([8]). (a) The mean exit times for f{x) = 
— X, d = l,a = 0.5, D = {—1, 1) and e = 0, 0.1, 0.5, 1. (b) Same as (a) except 
a = 1.5. 

in (j2ip become more "discontinuous" at the boundary as a — )• 0+ (more 
precisely, the derivative of u{x) goes to infinity faster near the boundary for 
smaller values of a). Our numerical results show that adding 0-U potential 
would cause the mean exit times be discontinuous at the boundary of the 
domain for all values of a in (0, 1) and any domain size. 

For large values of q (1 < a < 2), such as a = 1.5 as shown in Fig. [9|^a), 
the mean exit time shapes are similar to the parabolic shape as in the pure 
Gaussian noise case. Clearly, as the amount of Gaussian noises d increases 
keeping other factors fixed, the mean exit times decreases. As shown in 
the figure, the mean exit time is continuous at the boundary even in the 
absence of Gaussian noise {d = 0). From the numerical simulations (not 
shown here), we also find that the mean exit times are continuous at the 
boundary for 1 < a < 2. 

Next, we look at the effect of non-Gaussian noises by changing the pa- 
rameter £ while keeping /(x) = —x, d = 1 and D = (—1,1). It is obvious 
from the numerical results shown in Fig. [10] that the mean exit times de- 
creases when the amount of non-Gaussian noises e increases for all values 
of a in (0,2]. Due to the presence of significant Gaussian noises (d = 1), 
the shapes of the mean exit times u{x) are parabola-like for all parameter 
values shown in the figure. Keeping other parameters fixed, the effect of 
non-Gaussian noises on mean exit time is stronger when a is larger. It is 
consistent with the previous result shown in Fig. [7^a) that, for small do- 
mains, the mean exit times decreases as a increases. 
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Figure 11: The escape probability of the Levy motion, defined by the gener- 
ator in ([8|), out of the domain D landing in E. (a) The escape probabilities 
Pe{x) for the symmetric a-stable processes, i.e., d = 0, / = and e = 1, 
with a = 0.5,1,1.5,2, D = (—1,1) and E = [l,oo), and the analytical re- 
sults are shown by fine-dotted lines for all a values; (b) Same as (a) except 
D = (—4,4) and E = [4, oo); (c) Same as (a) except for the Levy motion 
with 0-U potential f{x) = —x, i.e., d = 0, f{x) = —x and e = 1; (d) Same 
as (c) except D = (—4,4) and E = [4, oo); (e) Same as (c) except with 
added Gaussian noise, i.e., d = 0.1, /(x) = —x and e = 1; (f) Same as (e) 
except D = (-4,4) and E = [4,oo). 



5.4 Escape probability 



In this section, we simulate the escape probability described by ([9]). In 
particular, for the special case of D = (a, b) and E = [b, oo), Eq. ([9]) becomes 



rl s^C 

>U^) + f{x)Pi,{x)-'-^ 
2 a 



1 



{x — a) 

Pe{x + y)- Pe{x) - /{|j^|<5}2/P^(x) 



+ 



1 



|i+« 



{b - xY 



Pe{x 
1 



a {b — xY 



,(28) 



for x G (a, b). The conditions for the escape probability outside the domain 
are Pe{x) = for a; G (— oo,a] and Pe{x) = 1 for x E [6, oo). 

First, we verify our numerical schemes by comparing with the analytical 
result of the escape probability for the symmetric Q-stable case (/ = 0, 
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d = 0,e = l) with D = (-1, 1) and E = [1, oo) [5] 



= HFT^ rib'-y')^-'dy, X e (-6,6). (29) 
[r(a/2)j^ 

As shown in Fig. [TlTa) and (b), the numerical results match with the 
analytical results given in Eq. (j29p . Due to the symmetry of the process and 
the domains, the escape probability takes the value of one-half when the 
starting point is the position of symmetry a: = 0. The escape probability is 
symmetric with respect to the point (0, 1/2), i.e., 

Pe{x) + Pe{-x) = 1. 

These remain true even we add Brownian noise and the 0-U potential to 
the process. 

Due to the symmetry, in the following discussion we focus on positive 
starting points in the domain, i.e., a; > 0. Figure [TT] shows that the proba- 
bility for the process to escape to the right of the domain is smaller when the 
value of a decreases. For a fixed positive x, the escape probability Pe{x) is 
the largest in for the case of Gaussian noise only (a = 2), keeping other fac- 
tors the same. This property is independent of the domain size or whether 
there exists a deterministic driving mechanism /. For larger domain sizes, 
as shown in Fig. [TTIfb.d.f). the escape probability tends to the value of equal 
chance 1/2 especially for small values of a. By comparing Fig. [TTT c) with 
Fig lll( a) or comparing Fig. [TTT d) with Fig. [TTIfb). we find that the effect 
of 0-U potential is reducing the escape probability for the same starting 
point X. The escape probability Pe{x) for the Brownian noise (q = 2) is no 
longer a straight line in the presence of 0-U potential. Again, for < a < 1, 
we find that the escape probability is discontinuous at the boundary of the 
domain when the SDE is driven by the 0-U potential and "pure" a-stable 
symmetric process, as demonstrated by the graphs of a = 0.5 in Fig. [TlTc) 
and (d). Adding Gaussian noises {d = 0.1) to the processes increases the 
chances of escape to the right, as shown in Fig. [TTIfe) and (f). 

Next, we consider SDE ([1]) driven by the double-well potential f{x) = x— 
x^. The corresponding deterministic dynamical system dXt = {Xt — Xf) dt 
has two stable states located at x = ±1 while is an unstable steady state. 
The double-well potential is well-known and widely used in phase transition 
studies. 

We investigate the likelihood of the stochastic process that starts within 
a bounded domain D C (— oo,0) and escapes and lands in the right-half 
line E = [0, oo) compared with that of escaping to the left of the bounded 
domain. We consider the bounded domains D that includes the left stable 
point x = —1 and the escape-target domain E containing the other stable 
point X = 1. When a process lands in E, in absence of the noises, it will 
be driven to and stay at the stable point x = 1. In other words, we try 
to examine the effect of noises on the likelihood of the transition from one 



21 




Figure 12: The escape probability of the Levy motion, defined by the gen- 
erator in ([8]), out of a bounded domain D landing in E, driven by the 
double- well potential f{x) = x — x"^ and noises, (a) The escape probability 
Pe{x) for d = 0,e = 1 D = (-1.1,0), E = [0,oo), a = 0.5,1,1.5,2; (b) 
Same as in (a) except D = (—2,0); (c) same as (a) except d = 0.1; (d) same 
as (b) except c? = 0.1. 
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stable state to the other. Figure [T2ra) shows the escape probabihty for 
D = (—1.1, 0), E = [0, oo) and a = 0.5, 1, 1.5, 2, when the stochastic effects 
are given by a-stable symmetric processes only (e = 1, and d = 0). The 
escape probability Pe{x) deviates more from a straight line as a decreases 
and it is smaller for smaller a when starting from x > 0.5. On the contrast, 
the probability is larger for smaller a when the starting point is close to the 
left boundary of the bounded domain. For the bigger domain D = (—2,0) 
shown in Figll2fb). the likelihood of escape to the right is more than a half 
for most of the starting points x > —1.5 and 1 < a < 2; the probability 
stays close to a half for most of the starting points when a = 0.5. Note that, 
in the case of a < 1, the probability is discontinuous at the left boundary 
of the domain. As shown in Fig]12fc) and (d), the differences in escape 
probabilities among different values of a become smaller when an amount 
of Gaussian noises is added {d = 0.1), but otherwise probabilities have the 
similar values and properties compared with those of d = 0. 

6 Conclusion 

In summary, we have developed an accurate numerical scheme for solving 
the mean first exit time and escape probability for SDKs with non-Gaussian 
Levy motions. We have analyzed the numerical error due to the singular 
nature of the Levy measure corresponding to jumps and accordingly, added 
a correction term to the numerical scheme. We have validated the numerical 
method by comparing with analytical and asymptotic solutions. For arbi- 
trary deterministic driving force, we have also given asymptotic solutions 
of the mean exit time when the pure jump measure in the Levy motion is 
small. 

Using both analytical and numerical results, we find that the mean exit 
time depends strongly on the domain size and the value of a in the a— stable 
Levy jump measure. For example, for a-stable Levy motion, the mean exit 
time can help us decide which of the two competing factors in a-stable Levy 
motion, the jump frequency or the jump size, is dominant. The small jumps 
with high frequency, corresponding to Levy motion with a closer to 2, make 
it easier to exit the small domains. On the other hand, it is easier to exit 
large domains for a-stable Levy motion with a closer to 0, which has the 
characteristics of the large jumps with low frequency. Another observation 
from the results is that, for smaller values of a (0 < a < 1), the mean exit 
time profiles are flat away from the boundary of the domain. This implies 
that the jump sizes of the processes are usually larger than the domain 
sizes, thus the mean exit times have small variations for different starting 
positions. 

The probability for the process to escape to the right of the domain is 
smaller when the value of a decreases. For a fixed positive x, the escape 
probability Pe{x) is the largest in for the case of Gaussian noise only (a = 2), 
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keeping other factors the same. This property is independent of the domain 
size or whether there exists a deterministic driving mechanism /. The escape 
probability is shown to vary significantly with the underlying vector field. 

The mean exit time and escape probability could become discontinu- 
ous at the boundary of the domain, when the process is subject to certain 
deterministic potential and the value of a is in (0, 1). 
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